

use final_data,clear
egen id=group(province)
duplicates drop  id year ,force
keep if year<2018
drop if co2==.
xtset id year
sum l_co2  p_new l_total l_forest l_population l_road l_area unemploy cpi gdp_a  phone  internet_use second_product l_patent pm l_fire 
des l_co2  p_new l_total l_forest l_population l_road l_area unemploy cpi gdp_a  phone  internet_usesecond_product l_patent  pm l_fire 


//unitroot test
xtunitroot fisher l_co2 ,pperron lags(2)
xtunitroot fisher p_new ,pperron lags(2)


collin p_new    l_total  l_forest l_population l_road l_area unemploy cpi gdp_a  phone  internet_use second_product  
//correlation

 pwcorr_a l_co2 pm p_new l_fire l_patent second_product l_forest l_population ///
 l_road l_area unemploy cpi gdp_a,star1(0.05) star5(0.01) star10(0.01)

 gen p_2=p_new^2
 //basic regression
 reg l_co2  p_new l_total  l_forest l_population l_road l_area  unemploy cpi gdp_a   
 est store reg1
 xtgls l_co2  p_new l_total  l_forest l_population l_road l_area  unemploy cpi  gdp_a i.year , panels(hetero) corr(ar1)
 est store reg2
  xtreg l_co2  p_new l_total  l_forest l_population l_road  l_area unemploy cpi  gdp_a ,fe
 est store reg3
 xtreg l_co2  p_new l_total  l_forest l_population l_road l_area  unemploy cpi   gdp_a ,re
 est store reg4
 xtreg l_co2  p_new l_total  l_forest l_population l_road  l_area unemploy cpi  gdp_a i.year ,fe
 est store reg5
hausman reg3 reg4
esttab reg1 reg2 reg3 reg4  reg5 , cells(b(star fmt(4)) t(par fmt(2))) /// 
legend label  varwidth(50) varlabels(_cons constant) mtitles("OLS" "FGLS" "FE" "RE" "FE" ) /// 
stats(N r2, fmt(0 3) label(Observation R-squared ))  star(* .10 ** .05 *** .01) /*table 4*/

//Lag effect
xtreg L(0/1).l_co2  p_new  l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a i.year,fe
 est store reg0
 xtreg L(0/1).l_co2  L(1).p_new  l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a i.year,fe
 est store reg1
 xtreg l_co2  L(0/1).p_new  l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a i.year,fe
 est store reg2
 xtreg l_co2 p_new D.p_new  l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a ,fe
est store reg3
 xtreg l_co2  p_new p_2 l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a i.year,fe
 est store reg4
  xtreg L(0/1).l_co2  p_new p_2 l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a i.year,fe
 est store reg5
esttab reg0 reg1 reg2 reg3   reg4 reg5, cells(b(star fmt(4)) t(par fmt(2))) /// 
legend label  varwidth(50) varlabels(_cons constant)  /// 
keep (L.l_co2 p_new p_2 L.p_new _cons D.p_new) ///
stats(N r2, fmt(0 3) label(Observation R-squared ))  star(* .10 ** .05 *** .01) /*table 4*/




//Dynamic model

tabstat l_patent,st(mean p50)
gen level=.
replace level=0 if l_patent<=8.731982
replace level=1 if l_patent>8.731982
gen l_patent_new=level*p_new
gen l_patent_new1=l_patent*p_new

xtreg l_co2  p_new   l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a i.year if level==0,fe
est store reg1
xtreg l_co2  p_new   l_total  l_forest l_population l_road l_area  unemploy cpi gdp_a i.year if level==1,fe
est store reg2

xtreg l_co2  p_new l_patent_new level l_total  l_forest l_population l_road l_area unemploy cpi gdp_a i.year ,fe
est store reg3
xtreg L(0/1).l_co2  p_new   l_total  l_forest l_population l_road l_area unemploy cpi  gdp_a i.year if level==0,fe
est store reg4
xtreg L(0/1).l_co2  p_new   l_total  l_forest l_population l_road l_area  unemploy cpi gdp_a i.year if level==1,fe
est store reg5

xtreg L(0/1).l_co2  p_new l_patent_new level l_total  l_forest l_population l_road l_area unemploy cpi gdp_a i.year ,fe
est store reg6
esttab reg2 reg1 reg3 reg5 reg4    reg6, cells(b(star fmt(4)) t(par fmt(2))) /// 
keep (L.l_co2 p_new  l_patent_new  level  _cons ) ///
legend label  varwidth(50) varlabels(_cons constant)  /// 
stats(N r2, fmt(0 4) label(Observation R-squared ))  star(* .10 ** .05 *** .01) /*table 4*/

pwcorr_a p_new pm l_fire ,star1(0.05) star5(0.01) star10(0.01)



//mediation effect
tabstat l_patent,st(mean p50)
gen level=.
replace level=0 if l_patent<=8.731982
replace level=1 if l_patent>8.731982
gen l_patent_new=level*p_new
xtreg l_co2  p_new    l_total  l_forest l_population l_road l_area    gdp_a i.year  if level==0,fe
est store reg1
xtreg l_co2  p_new    l_total  l_forest l_population l_road l_area    gdp_a i.year if level==1,fe
est store reg2
xtreg l_co2  p_new l_patent_new  level l_total  l_forest l_population l_road l_area    gdp_a i.year ,fe
est store reg3
xtreg L(0/1).l_co2  p_new    l_total  l_forest l_population l_road l_area    gdp_a i.year  if level==0,fe
est store reg4
xtreg L(0/1).l_co2  p_new    l_total  l_forest l_population l_road l_area    gdp_a i.year if level==1,fe
est store reg5
xtreg L(0/1).l_co2  p_new l_patent_new  level l_total  l_forest l_population l_road l_area    gdp_a i.year ,fe
est store reg6
esttab reg1 reg2  reg3  reg4 reg5 reg6, cells(b(star fmt(4)) t(par fmt(2))) /// 
legend label  varwidth(50) varlabels(_cons constant)  /// 
stats(N r2, fmt(0 3) label(Observation R-squared ))  star(* .10 ** .05 *** .01) /*table 4*/

//figure of mediation
xtreg l_co2  c.p_new##c.level    l_total  l_forest l_population l_road l_area    gdp_a i.year , fe
est sto regression
foreach v of var p_new level {
  su `v' if e(sample)
  local low_`v'=r(mean)-r(sd)
  local high_`v'=r(mean)+r(sd)
}

est restore regression //
margins,at(p_new=(`low_p_new' `high_p_new') ///
         level = (`low_level' `high_level'))  //
marginsplot , xlabel(`low_p_new' "Low IV" `high_p_new' "High IV")  ///
              ytitle("Dependent variable")       ///
              ylabel(, angle(horizontal) nogrid) ///
              legend(position(7) col(1) stack)   ///
              title("") noci

///IV
xtivreg l_co2  (p_new=pm l_fire ) l_total  l_forest l_population l_road l_area  unemploy cpi  gdp_a   ,fe 
est store reg3
xtoverid
dmexogxt
xtivreg L(0/1).l_co2  (p_new=pm l_fire ) l_total  l_forest l_population l_road l_area  unemploy cpi  gdp_a   ,fe 
est store reg4
xtoverid
dmexogxt
xtreg l_co2  p_new    l_total  l_forest l_population l_road l_area unemploy cpi gdp_a phone  internet_use second_product i.year ,fe
est store reg1
reg l_co2  p_new    l_total  l_forest l_population l_road l_area unemploy cpi gdp_a phone  internet_use second_product
est store reg2
esttab reg3 reg4   reg2 reg1 , cells(b(star fmt(4)) t(par fmt(2))) /// 
legend label  varwidth(50) varlabels(_cons constant) mtitles("全样本" "2010" "2012" "2014" "2016" "FE") /// 
stats(N r2, fmt(0 3) label(Observation R-squared ))  star(* .10 ** .05 *** .01) /*table 4*/



